home *** CD-ROM | disk | FTP | other *** search
/ The Atari Compendium / The Atari Compendium (Toad Computers) (1994).iso / files / prgtools / gnustuff / tos / pmathlib / pmlsrc18.zoo / pmlsrc / sqrt.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-03-27  |  9.1 KB  |  426 lines
  1. /************************************************************************
  2.  *                                    *
  3.  *                N O T I C E                *
  4.  *                                    *
  5.  *            Copyright Abandoned, 1987, Fred Fish        *
  6.  *                                    *
  7.  *    This previously copyrighted work has been placed into the    *
  8.  *    public domain by the author (Fred Fish) and may be freely used    *
  9.  *    for any purpose, private or commercial.  I would appreciate    *
  10.  *    it, as a courtesy, if this notice is left in all copies and    *
  11.  *    derivative works.  Thank you, and enjoy...            *
  12.  *                                    *
  13.  *    The author makes no warranty of any kind with respect to this    *
  14.  *    product and explicitly disclaims any implied warranties of    *
  15.  *    merchantability or fitness for any particular purpose.        *
  16.  *                                    *
  17.  ************************************************************************
  18.  */
  19.  
  20.  
  21. /*
  22.  *  FUNCTION
  23.  *
  24.  *    sqrt   double precision square root
  25.  *
  26.  *  KEY WORDS
  27.  *
  28.  *    sqrt
  29.  *    machine independent routines
  30.  *    math libraries
  31.  *
  32.  *  DESCRIPTION
  33.  *
  34.  *    Returns double precision square root of double precision
  35.  *    floating point argument.
  36.  *
  37.  *  USAGE
  38.  *
  39.  *    double sqrt (x)
  40.  *    double x;
  41.  *
  42.  *  REFERENCES
  43.  *
  44.  *    Fortran IV-PLUS user's guide, Digital Equipment Corp. pp B-7.
  45.  *
  46.  *    Computer Approximations, J.F. Hart et al, John Wiley & Sons,
  47.  *    1968, pp. 89-96.
  48.  *
  49.  *  RESTRICTIONS
  50.  *
  51.  *    The relative error is 10**(-30.1) after three applications
  52.  *    of Heron's iteration for the square root.
  53.  *
  54.  *    However, this assumes exact arithmetic in the iterations
  55.  *    and initial approximation.  Additional errors may occur
  56.  *    due to truncation, rounding, or machine precision limits.
  57.  *    
  58.  *  PROGRAMMER
  59.  *
  60.  *    Fred Fish
  61.  *
  62.  *  INTERNALS
  63.  *
  64.  *    Computes square root by:
  65.  *
  66.  *      (1)    Range reduction of argument to [0.5,1.0]
  67.  *        by application of identity:
  68.  *
  69.  *        sqrt(x)  =  2**(k/2) * sqrt(x * 2**(-k))
  70.  *
  71.  *        k is the exponent when x is written as
  72.  *        a mantissa times a power of 2 (m * 2**k).
  73.  *        It is assumed that the mantissa is
  74.  *        already normalized (0.5 =< m < 1.0).
  75.  *
  76.  *      (2)    An approximation to sqrt(m) is obtained
  77.  *        from:
  78.  *
  79.  *        u = sqrt(m) = (P0 + P1*m) / (Q0 + Q1*m)
  80.  *
  81.  *        P0 = 0.594604482
  82.  *        P1 = 2.54164041
  83.  *        Q0 = 2.13725758
  84.  *        Q1 = 1.0
  85.  *
  86.  *        (coefficients from HART table #350 pg 193)
  87.  *
  88.  *      (3)    Three applications of Heron's iteration are
  89.  *        performed using:
  90.  *
  91.  *        y[n+1] = 0.5 * (y[n] + (m/y[n]))
  92.  *
  93.  *        where y[0] = u = approx sqrt(m)
  94.  *
  95.  *      (4)    If the value of k was odd then y is either
  96.  *        multiplied by the square root of two or
  97.  *        divided by the square root of two for k positive
  98.  *        or negative respectively.  This rescales y
  99.  *        by multiplying by 2**frac(k/2).
  100.  *
  101.  *      (5)    Finally, y is rescaled by int(k/2) which
  102.  *        is equivalent to multiplication by 2**int(k/2).
  103.  *
  104.  *        The result of steps 4 and 5 is that the value
  105.  *        of y between 0.5 and 1.0 has been rescaled by
  106.  *        2**(k/2) which removes the original rescaling
  107.  *        done prior to finding the mantissa square root.
  108.  *
  109.  *  NOTES
  110.  *
  111.  *    The Convergent Technologies compiler optimizes division
  112.  *    by powers of two to "arithmetic shift right" instructions.
  113.  *    This is ok for positive integers but gives different
  114.  *    results than other C compilers for negative integers.
  115.  *    For example, (-1)/2 is -1 on the CT box, 0 on every other
  116.  *    machine I tried.
  117.  *
  118.  */
  119.  
  120. /*
  121.  *  MODIFICATIONS
  122.  *
  123.  *    This routine modified to use polynomial, instead of rational,
  124.  *    approximation with coefficients
  125.  *
  126.  *        P0  0.22906994529e+00
  127.  *        P1  0.13006690496e+01
  128.  *        P2 -0.90932104982e+00
  129.  *        P3  0.50104207633e+00
  130.  *        P4 -0.12146838249e+00
  131.  *
  132.  *    as given by Hart (op. cit.) in SQRT 0132.  This approximation
  133.  *    gives around 5 digits correct.  Two applications of Heron's
  134.  *    approximation will give more then practically achievable
  135.  *    13-15 decimal digits.  Multiplications by powers of 2 are
  136.  *    replaced by explicit calls to ldexp.
  137.  *
  138.  *    Michal Jaegermann, 24 October 1990
  139.  */
  140.  
  141. #if !defined (__M68881__) && !defined (sfp004)
  142.  
  143. #include <stdio.h>
  144. #include <math.h>
  145. #include "pml.h"
  146.  
  147. #ifdef OLD
  148. #define P0 0.594604482            /* Approximation coeff    */
  149. #define P1 2.54164041            /* Approximation coeff    */
  150. #define Q0 2.13725758            /* Approximation coeff    */
  151. #define Q1 1.0                /* Approximation coeff    */
  152.  
  153. #define ITERATIONS 3            /* Number of iterations    */
  154.  
  155. #endif
  156.  
  157. #define P0  0.22906994529e+00        /* Hart SQRT 0132 */
  158. #define P1  0.13006690496e+01
  159. #define P2 -0.90932104982e+00
  160. #define P3  0.50104207633e+00
  161. #define P4 -0.12146838249e+00
  162.  
  163. static char funcname[] = "sqrt";
  164.  
  165. double sqrt (x)
  166. double x;
  167. {
  168. #ifdef OLD
  169.     int k;
  170.     register int bugfix;
  171.     register int kmod2;
  172.     register int count;
  173.     int exponent;
  174.     double y;
  175. #else
  176.     int k;
  177. #endif
  178.     double m;
  179.     double u;
  180.     struct exception xcpt;
  181.  
  182.     if (x == 0.0) {
  183.     xcpt.retval = 0.0;
  184.     } else if (x < 0.0) {
  185.     xcpt.type = DOMAIN;
  186.     xcpt.name = funcname;
  187.     xcpt.arg1 = x;
  188.     if (!matherr (&xcpt)) {
  189.         fprintf (stderr, "%s: DOMAIN error\n", funcname);
  190.         errno = EDOM;
  191.         xcpt.retval = 0.0;
  192.     }
  193.     } else {
  194.     m = frexp (x, &k);
  195. #ifdef OLD
  196.     u = (P0 + (P1 * m)) / (Q0 + (Q1 * m));
  197.     for (count = 0, y = u; count < ITERATIONS; count++) {
  198.         y = 0.5 * (y + (m / y));
  199.     }
  200.     if ((kmod2 = (k % 2)) < 0) {
  201.         y /= SQRT2;
  202.     } else if (kmod2 > 0) {
  203.         y *= SQRT2;
  204.     }
  205.     bugfix = 2;
  206.     xcpt.retval = ldexp (y, k/bugfix);
  207. #else
  208.     u = (((P4 * m + P3) * m + P2) * m + P1) * m + P0;
  209.     u = ldexp((u + (m / u)), -1);    /* Heron's iteration */
  210.     u += m / u;            /* and a part of the second one */
  211.     if (k & 1) {
  212.         u *= SQRT2;
  213.     }
  214.     /*
  215.      * here we rely on the fact that -3/2 and (-3 >> 1)
  216.      * do give different results
  217.      */
  218.     xcpt.retval = ldexp (u, (k >> 1) - 1);
  219. #endif
  220.     }
  221.     return (xcpt.retval);
  222. }
  223.  
  224. double hypot(double x, double y)
  225. {
  226.     return sqrt(x*x + y*y);
  227. }
  228. #endif    /*  __M68881__, sfp004    */
  229.  
  230. #ifdef    sfp004
  231.  
  232. __asm("
  233.  
  234. comm =     -6
  235. resp =    -16
  236. zahl =      0
  237.  
  238. ");    /* end asm    */
  239.  
  240. #endif    sfp004
  241. #if defined (__M68881__) || defined (sfp004)
  242.  
  243.     __asm(".text; .even");
  244.  
  245. # ifdef    ERROR_CHECK
  246.  
  247.     __asm("
  248.  
  249. _Overflow:
  250.     .ascii \"OVERFLOW\\0\"
  251. _Domain:
  252.     .ascii \"DOMAIN\\0\"
  253. _Error_String:
  254.     .ascii \"sqrt: %s error\\n\\0\"
  255. .even
  256.  
  257. | m.ritzert 7.12.1991
  258. | ritzert@dfg.dbp.de
  259. |
  260. |    /* NAN  = {7fffffff,ffffffff}        */
  261. |    /* +Inf = {7ff00000,00000000}        */
  262. |    /* -Inf = {fff00000,00000000}        */
  263. |    /* MAX_D= {7fee42d1,30773b76}        */
  264. |    /* MIN_D= {ffee42d1,30773b76}        */
  265.  
  266. .even
  267. double_max:
  268.     .long    0x7fee42d1
  269.     .long    0x30273b76
  270. double_min:
  271.     .long    0xffee42d1
  272.     .long    0x30273b76
  273. NaN:
  274.     .long    0x7fffffff
  275.     .long    0xffffffff
  276. p_Inf:
  277.     .long    0x7ff00000
  278.     .long    0x00000000
  279. m_Inf:
  280.     .long    0xfff00000
  281.     .long    0x00000000
  282.     ");    /* end asm    */
  283. # endif    ERROR_CHECK
  284.  
  285.     __asm(".even
  286. .globl _hypot
  287. .globl _sqrt
  288. _sqrt:
  289.     ");    /* end asm    */
  290.  
  291. #endif    /* __M68881__ || sfp004    */
  292. #ifdef    __M68881__
  293.  
  294.     __asm("
  295.     fsqrtd    a7@(4), fp0    | sqrt
  296.     fmoved    fp0,a7@-    | push result
  297.     moveml    a7@+,d0-d1    | return_value
  298.     ");    /* end asm    */
  299.  
  300. #endif    __M68881__
  301. #ifdef    sfp004
  302.     __asm("
  303.     lea    0xfffa50,a0
  304.     movew    #0x5404,a0@(comm)    | specify function
  305.     cmpiw    #0x8900,a0@(resp)    | check
  306.     movel    a7@(4),a0@        | load arg_hi
  307.     movel    a7@(8),a0@        | load arg_low
  308.     movew    #0x7400,a0@(comm)    | result to d0
  309.     .long    0x0c688900, 0xfff067f8    | wait
  310.     movel    a0@,d0
  311.     movel    a0@,d1
  312.     ");    /* end asm    */
  313.  
  314. #endif    sfp004
  315. #if defined (__M68881__) || defined (sfp004)
  316. # ifdef    ERROR_CHECK
  317.     __asm("
  318. err:
  319.     lea    double_max,a0    |
  320.     swap    d0        | exponent into lower word
  321.     cmpw    a0@(16),d0    | == NaN ?
  322.     beq    error_nan    |
  323.     cmpw    a0@(24),d0    | == + Infinity ?
  324.     beq    error_plus    |
  325.     swap    d0        | result ok,
  326.     rts            | restore d0
  327. ");
  328. #ifndef    __MSHORT__
  329. __asm("
  330. error_plus:
  331.     swap    d0
  332.     moveml    d0-d1,a7@-
  333.     movel    #63,_errno    | NAN => errno = EDOM
  334.     pea    _Overflow    | for printf
  335.     bra    error_exit    |
  336. error_nan:
  337.     moveml    a0@(24),d0-d1    | result = +inf
  338.     moveml    d0-d1,a7@-
  339.     movel    #62,_errno    | NAN => errno = EDOM
  340.     pea    _Domain        | for printf
  341. ");
  342. #else    __MSHORT__
  343. __asm("
  344. error_plus:
  345.     swap    d0
  346.     moveml    d0-d1,a7@-
  347.     movew    #63,_errno    | NAN => errno = EDOM
  348.     pea    _Overflow    | for printf
  349.     bra    error_exit    |
  350. error_nan:
  351.     moveml    a0@(24),d0-d1    | result = +inf
  352.     moveml    d0-d1,a7@-
  353.     movew    #62,_errno    | NAN => errno = EDOM
  354.     pea    _Domain        | for printf
  355. ");
  356. #endif    __MSHORT__
  357. __asm("
  358. error_exit:
  359.     pea    _Error_String    |
  360.     pea    __iob+52    |
  361.     jbsr    _fprintf    |
  362.     addl    #12,a7        |
  363.     moveml    a7@+,d0-d1
  364. ");
  365. # endif    ERROR_CHECK
  366. __asm("
  367.     rts
  368.  
  369. .even
  370. _hypot:
  371.     ");
  372. #endif /* __M68881__ || sfp004    */
  373. #ifdef __M68881__
  374. __asm("
  375.     fmoved    a7@(4),fp0    |
  376.     fmulx    fp0,fp0        | x**2
  377.     fmoved    a7@(12),fp1    |
  378.     fmulx    fp1,fp1        | y**2
  379.     faddx    fp1,fp0        |
  380.     fsqrtx    fp0,fp0        | sqrt( x**2 + y**2 )
  381.     fmoved    fp0,a7@-    |
  382.     moveml    a7@+,d0-d1    | return arg
  383. ");
  384. #endif    __M68881__
  385. #ifdef    sfp004
  386. __asm("
  387.     lea    0xfffa50,a0
  388.  
  389.     movew    #0x5400,a0@(comm)    | load fp0
  390.     .long    0x0c688900, 0xfff067f8
  391.     movel    a7@(4),a0@        | load arg_hi
  392.     movel    a7@(8),a0@        | load arg_low
  393.  
  394.     movew    #0x5480,a0@(comm)    | load fp1
  395.     .long    0x0c688900, 0xfff067f8
  396.     movel    a7@(12),a0@        | load arg_hi
  397.     movel    a7@(16),a0@        | load arg_low
  398.  
  399.     movew    #0x0023,a0@(comm)
  400.     .word    0x4a68,0xfff0,0x6bfa    | test
  401.  
  402.     movew    #0x04a3,a0@(comm)
  403.     .word    0x4a68,0xfff0,0x6bfa    | test
  404.  
  405.     movew    #0x0422,a0@(comm)    | fp0 = fp0 + fp1    
  406.     .word    0x4a68,0xfff0,0x6bfa    | test
  407.  
  408.     movew    #0x0004,a0@(comm)    | sqrt(fp0)
  409.     .word    0x4a68,0xfff0,0x6bfa    | test
  410.  
  411.     movew    #0x7400,a0@(comm)    | result to d0/d1
  412.     .long    0x0c688900, 0xfff067f8
  413.     movel    a0@(zahl),d0
  414.     movel    a0@(zahl),d1
  415. ");
  416. #endif    sfp004
  417. #if ( defined (__M68881__) || defined (sfp004) ) && defined (ERROR_CHECK)
  418.  
  419. __asm("bra err");
  420.  
  421. #else    ERROR_CHECK
  422.  
  423. __asm("rts");
  424.  
  425. #endif    ERROR_CHECK
  426.